#!/bin/csh -f
#
# qsumm: Script to prepare summary file from Aqua output
# =============================================================
# Ton Rullmann, Mar 3, 1997 - Aqua 1.2
# =============================================================
# The variable 'aquaroot' must exist and point to the Aqua top directory.
# The 'aqsetup' script should reside in that directory.


if ( $#argv < 2 ) then
    echo "usage: qsumm file1 file2 [numlin]"
    echo "    file1: violation data"
    echo "    file2: information content data"
    echo "   numlin: maximum number of lines shown in sorted list of largest"
    echo "           violations"
    exit( 1 )
endif
set numlin = 3
if ( $#argv >= 3 ) set numlin = $3

# flag for removing empty chain ID's - can't be used for files in which 
# non-empty chain ID's occur 
set oldstyle = 0
if ( $#argv == 4 ) set oldstyle = 1

source $aquaroot/aqsetup

# Remove temporary files
if (-e aqsumm.outu) \rm aqsumm.outu
if (-e aqsumm.outl) \rm aqsumm.outl

echo "SUMMARY  OF  AQUA  RESULTS  "`date`  "(Aqua version $AQUAVERS)"
gawk 'NR==1{print "DATA: ",$(NF-1),$NF; exit}' $2
echo "from files:"
echo '<'$1'> and <'$2'>'

#################################
# Restraints per residue and type
#################################

echo " "
echo "Information content"
echo "==================="

# Extract the number of structures
set nresid = `qdbext -v "b=general" -v "i=nresidues" -v "hdr=none" $2`
set nrestr = `qdbext -v "b=general" -v "i=nrestraints" -v "hdr=none" $2`
qdbext -v "b=NUMBER_OF_RESTRAINTS_per_RESIDUE" -v "hdr=none" $2 | \
    gawk -f $AQUASCRIPTSDIR/qsumm_aux1 -v "nresid=$nresid" -v "nrestr=$nrestr"
qdbext -v "b=NUMBER_OF_RESTRAINTS_per_RANGE_TYPE" -v "hdr=none" $2 | gawk \
	' BEGIN { print "\nNumber of restraints per range type:" } \
	{ printf( "%12s %5i\n", $1, $NF ) }'

##########################
# Extract violation data 
##########################

# Extract the number of structures
set nstruc = `qdbext -v "b=general" -v "i=nstructures" -v "hdr=none" $1`

# Extract the data record with all restraint violations
qdbext -v "b=RESTRAINT_VIOLATIONS" -v "hdr=none" $1 > aqsumm.all

# Extract the data record with restraint violation summary per structure
qdbext -v "b=RESTRAINT_VIOLATION_SUMMARY_per_STRUCTURE" -v "hdr=none" $1 > aqsumm.str

##########################
# Loop over all structures -> print largest and average violations
##########################

set istruc = 1
set off = 10
while ( $istruc <= $nstruc )
    echo " "
    echo "Model $istruc "
    echo "========"
    @ i = $istruc + $off
    @ j = $i + 1

# Upper bound violations
# Sort on the corresponding column of violation data, print the largest ones (print the average)
    sort -r -n --key=$j,$j aqsumm.all | \
        gawk -f $AQUASCRIPTSDIR/qsumm_aux2 -v "num=$j" -v "aver=0" -v "numlin=$numlin" -v "oldstyle=$oldstyle"
    echo "Largest UPPER BOUND violations:"
    if ( -e aqsumm.tmp ) then
	cat aqsumm.tmp*
# Add the list of largest violations to aqsumm.outu
	gawk -v "struc=$istruc" \
	    '{print $0,"(model",struc")"}' \
	    aqsumm.tmp >> aqsumm.outu 
	\rm aqsumm.tmp*
    endif

# Lower bound violations
# Sort on the corresponding column of violation data, print the largest ones (print the average)
    sort -n --key=$j,$j aqsumm.all |\
        gawk -f $AQUASCRIPTSDIR/qsumm_aux2 -v "num=$j" -v "aver=0" -v "numlin=$numlin" -v "oldstyle=$oldstyle"
    echo "Largest LOWER BOUND violations:"
    if ( -e aqsumm.tmp ) then
	cat aqsumm.tmp*
# Add the list of largest violations to aqsumm.outl
	gawk -v "struc=$istruc" \
	    '{print $0,"(model",struc")"}' \
	    aqsumm.tmp >> aqsumm.outl 
	\rm aqsumm.tmp*
    endif

# Print the data from the SUMMARY_PER_STRUCTURE for this structure
    gawk -f $AQUASCRIPTSDIR/qsumm_aux3 -v "sel=$istruc" aqsumm.str 
    @ istruc++
end

##########################
# Overall violation data -> print largest and average violations
##########################

# Sort the merge files with largest violation data and print the largest

set i = 8
set j = 9
if ( $oldstyle ) then
    set i = 6
    set j = 7
endif

echo "Overall"
echo "======="
echo "Largest UPPER BOUND violations:"
if ( -e aqsumm.outu ) then
    sort -r -n --key=$j,$j aqsumm.outu | \
        gawk -v "num=$j" -v "numlin=$numlin"\
	'NR==numlin { exitval = $num } \
	$num < exitval { exit } \
	{ print $0 }' 
endif
echo "Largest LOWER BOUND violations:"
if ( -e aqsumm.outl ) then
    sort -n --key=$j,$j aqsumm.outl | \
        gawk -v "num=$j" -v "numlin=$numlin"\
    	'NR==numlin { exitval = $num } \
	$num > exitval { exit } \
	{ print $0 }' 
endif

# Print the data from the OVERALL SUMMARY
qdbext -v "b=RESTRAINT_VIOLATION_SUMMARY_OVERALL" -v "hdr=none" $1 | \
    gawk -f $AQUASCRIPTSDIR/qsumm_aux3 -v sel=1

##########################
# Per residue range violation data -> print largest and average violations
##########################

echo "Per residue range"
echo "================="
qdbext -v "b=RESTRAINT_VIOLATION_SUMMARY_per_RANGE_TYPE_for_ALL_STRUCTURES" -v "hdr=none" $1 | \
    gawk -f $AQUASCRIPTSDIR/qsumm_aux3 -v sel=0

##########################
# Finish
##########################

echo "For more details:"
echo "================="
echo "file <"$1">: detailed violation data"
echo "     (e.g. per residue, per restraint, per range type)"
echo "file <"$2">: detailed information content data"
echo "     (e.g. per residue, per range type)"
echo "These files can be parsed using the <qdbext> command."
echo " "

# Remove temporary files
#\rm aqsumm.*
